home *** CD-ROM | disk | FTP | other *** search
/ IRIX 6.2 Development Libraries / SGI IRIX 6.2 Development Libraries.iso / dist / complib.idb / usr / share / catman / p_man / cat3 / complib / dposvx.z / dposvx
Text File  |  1996-03-14  |  10KB  |  265 lines
  1.  
  2.  
  3.  
  4. DDDDPPPPOOOOSSSSVVVVXXXX((((3333FFFF))))                                                          DDDDPPPPOOOOSSSSVVVVXXXX((((3333FFFF))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      DPOSVX - use the Cholesky factorization A = U**T*U or A = L*L**T to
  10.      compute the solution to a real system of linear equations  A * X = B,
  11.  
  12. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  13.      SUBROUTINE DPOSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, EQUED, S, B,
  14.                         LDB, X, LDX, RCOND, FERR, BERR, WORK, IWORK, INFO )
  15.  
  16.          CHARACTER      EQUED, FACT, UPLO
  17.  
  18.          INTEGER        INFO, LDA, LDAF, LDB, LDX, N, NRHS
  19.  
  20.          DOUBLE         PRECISION RCOND
  21.  
  22.          INTEGER        IWORK( * )
  23.  
  24.          DOUBLE         PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
  25.                         BERR( * ), FERR( * ), S( * ), WORK( * ), X( LDX, * )
  26.  
  27. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  28.      DPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
  29.      compute the solution to a real system of linear equations
  30.         A * X = B, where A is an N-by-N symmetric positive definite matrix and
  31.      X and B are N-by-NRHS matrices.
  32.  
  33.      Error bounds on the solution and a condition estimate are also provided.
  34.  
  35.  
  36. DDDDEEEESSSSCCCCRRRRIIIIPPPPTTTTIIIIOOOONNNN
  37.      The following steps are performed:
  38.  
  39.      1. If FACT = 'E', real scaling factors are computed to equilibrate
  40.         the system:
  41.            diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
  42.         Whether or not the system will be equilibrated depends on the
  43.         scaling of the matrix A, but if equilibration is used, A is
  44.         overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
  45.  
  46.      2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
  47.         factor the matrix A (after equilibration if FACT = 'E') as
  48.            A = U**T* U,  if UPLO = 'U', or
  49.            A = L * L**T,  if UPLO = 'L',
  50.         where U is an upper triangular matrix and L is a lower triangular
  51.         matrix.
  52.  
  53.      3. The factored form of A is used to estimate the condition number
  54.         of the matrix A.  If the reciprocal of the condition number is
  55.         less than machine precision, steps 4-6 are skipped.
  56.  
  57.      4. The system of equations is solved for X using the factored form
  58.         of A.
  59.  
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. DDDDPPPPOOOOSSSSVVVVXXXX((((3333FFFF))))                                                          DDDDPPPPOOOOSSSSVVVVXXXX((((3333FFFF))))
  71.  
  72.  
  73.  
  74.      5. Iterative refinement is applied to improve the computed solution
  75.         matrix and calculate error bounds and backward error estimates
  76.         for it.
  77.  
  78.      6. If equilibration was used, the matrix X is premultiplied by
  79.         diag(S) so that it solves the original system before
  80.         equilibration.
  81.  
  82.  
  83. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  84.      FACT    (input) CHARACTER*1
  85.              Specifies whether or not the factored form of the matrix A is
  86.              supplied on entry, and if not, whether the matrix A should be
  87.              equilibrated before it is factored.  = 'F':  On entry, AF
  88.              contains the factored form of A.  If EQUED = 'Y', the matrix A
  89.              has been equilibrated with scaling factors given by S.  A and AF
  90.              will not be modified.  = 'N':  The matrix A will be copied to AF
  91.              and factored.
  92.              = 'E':  The matrix A will be equilibrated if necessary, then
  93.              copied to AF and factored.
  94.  
  95.      UPLO    (input) CHARACTER*1
  96.              = 'U':  Upper triangle of A is stored;
  97.              = 'L':  Lower triangle of A is stored.
  98.  
  99.      N       (input) INTEGER
  100.              The number of linear equations, i.e., the order of the matrix A.
  101.              N >= 0.
  102.  
  103.      NRHS    (input) INTEGER
  104.              The number of right hand sides, i.e., the number of columns of
  105.              the matrices B and X.  NRHS >= 0.
  106.  
  107.      A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  108.              On entry, the symmetric matrix A, except if FACT = 'F' and EQUED
  109.              = 'Y', then A must contain the equilibrated matrix
  110.              diag(S)*A*diag(S).  If UPLO = 'U', the leading N-by-N upper
  111.              triangular part of A contains the upper triangular part of the
  112.              matrix A, and the strictly lower triangular part of A is not
  113.              referenced.  If UPLO = 'L', the leading N-by-N lower triangular
  114.              part of A contains the lower triangular part of the matrix A, and
  115.              the strictly upper triangular part of A is not referenced.  A is
  116.              not modified if FACT = 'F' or 'N', or if FACT = 'E' and EQUED =
  117.              'N' on exit.
  118.  
  119.              On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
  120.              diag(S)*A*diag(S).
  121.  
  122.      LDA     (input) INTEGER
  123.              The leading dimension of the array A.  LDA >= max(1,N).
  124.  
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. DDDDPPPPOOOOSSSSVVVVXXXX((((3333FFFF))))                                                          DDDDPPPPOOOOSSSSVVVVXXXX((((3333FFFF))))
  137.  
  138.  
  139.  
  140.      AF      (input or output) DOUBLE PRECISION array, dimension (LDAF,N)
  141.              If FACT = 'F', then AF is an input argument and on entry contains
  142.              the triangular factor U or L from the Cholesky factorization A =
  143.              U**T*U or A = L*L**T, in the same storage format as A.  If EQUED
  144.              .ne. 'N', then AF is the factored form of the equilibrated matrix
  145.              diag(S)*A*diag(S).
  146.  
  147.              If FACT = 'N', then AF is an output argument and on exit returns
  148.              the triangular factor U or L from the Cholesky factorization A =
  149.              U**T*U or A = L*L**T of the original matrix A.
  150.  
  151.              If FACT = 'E', then AF is an output argument and on exit returns
  152.              the triangular factor U or L from the Cholesky factorization A =
  153.              U**T*U or A = L*L**T of the equilibrated matrix A (see the
  154.              description of A for the form of the equilibrated matrix).
  155.  
  156.      LDAF    (input) INTEGER
  157.              The leading dimension of the array AF.  LDAF >= max(1,N).
  158.  
  159.      EQUED   (input or output) CHARACTER*1
  160.              Specifies the form of equilibration that was done.  = 'N':  No
  161.              equilibration (always true if FACT = 'N').
  162.              = 'Y':  Equilibration was done, i.e., A has been replaced by
  163.              diag(S) * A * diag(S).  EQUED is an input argument if FACT = 'F';
  164.              otherwise, it is an output argument.
  165.  
  166.      S       (input or output) DOUBLE PRECISION array, dimension (N)
  167.              The scale factors for A; not accessed if EQUED = 'N'.  S is an
  168.              input argument if FACT = 'F'; otherwise, S is an output argument.
  169.              If FACT = 'F' and EQUED = 'Y', each element of S must be
  170.              positive.
  171.  
  172.      B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
  173.              On entry, the N-by-NRHS right hand side matrix B.  On exit, if
  174.              EQUED = 'N', B is not modified; if EQUED = 'Y', B is overwritten
  175.              by diag(S) * B.
  176.  
  177.      LDB     (input) INTEGER
  178.              The leading dimension of the array B.  LDB >= max(1,N).
  179.  
  180.      X       (output) DOUBLE PRECISION array, dimension (LDX,NRHS)
  181.              If INFO = 0, the N-by-NRHS solution matrix X to the original
  182.              system of equations.  Note that if EQUED = 'Y', A and B are
  183.              modified on exit, and the solution to the equilibrated system is
  184.              inv(diag(S))*X.
  185.  
  186.      LDX     (input) INTEGER
  187.              The leading dimension of the array X.  LDX >= max(1,N).
  188.  
  189.      RCOND   (output) DOUBLE PRECISION
  190.              The estimate of the reciprocal condition number of the matrix A
  191.              after equilibration (if done).  If RCOND is less than the machine
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. DDDDPPPPOOOOSSSSVVVVXXXX((((3333FFFF))))                                                          DDDDPPPPOOOOSSSSVVVVXXXX((((3333FFFF))))
  203.  
  204.  
  205.  
  206.              precision (in particular, if RCOND = 0), the matrix is singular
  207.              to working precision.  This condition is indicated by a return
  208.              code of INFO > 0, and the solution and error bounds are not
  209.              computed.
  210.  
  211.      FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
  212.              The estimated forward error bound for each solution vector X(j)
  213.              (the j-th column of the solution matrix X).  If XTRUE is the true
  214.              solution corresponding to X(j), FERR(j) is an estimated upper
  215.              bound for the magnitude of the largest element in (X(j) - XTRUE)
  216.              divided by the magnitude of the largest element in X(j).  The
  217.              estimate is as reliable as the estimate for RCOND, and is almost
  218.              always a slight overestimate of the true error.
  219.  
  220.      BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
  221.              The componentwise relative backward error of each solution vector
  222.              X(j) (i.e., the smallest relative change in any element of A or B
  223.              that makes X(j) an exact solution).
  224.  
  225.      WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
  226.  
  227.      IWORK   (workspace) INTEGER array, dimension (N)
  228.  
  229.      INFO    (output) INTEGER
  230.              = 0: successful exit
  231.              < 0: if INFO = -i, the i-th argument had an illegal value
  232.              > 0: if INFO = i, and i is
  233.              <= N: the leading minor of order i of A is not positive definite,
  234.              so the factorization could not be completed, and the solution and
  235.              error bounds could not be computed.  = N+1: RCOND is less than
  236.              machine precision.  The factorization has been completed, but the
  237.              matrix is singular to working precision, and the solution and
  238.              error bounds have not been computed.
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.  
  248.  
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.